clear all;
clc;

% ----------------------- Importing Data ----------------------------------

patnum = 5; % 1-5
th_val = -10; % 0, -10, and -20

cordgeom = load('patient_cordgeom.txt');

epi_loctag = ['_d',num2str(cordgeom(patnum,8)),'_t',num2str(th_val)];
epi_loctag = strrep(epi_loctag,'.','p');
epi_loctag = strrep(epi_loctag,'-','n');

sub_loctag = ['_d',num2str(cordgeom(patnum,9)),'_t',num2str(th_val)];
sub_loctag = strrep(sub_loctag,'.','p');
sub_loctag = strrep(sub_loctag,'-','n');

cd(['patient',num2str(patnum),'_thresholds']);
noaxons = 200;
epi_thdata = load(['neg_epi_','bp_dc_act_p',num2str(patnum),epi_loctag,'.txt']);
epi_thdata_Lbox = zeros(10,10);
epi_thdata_Rbox = zeros(10,10);
epi_drthresh = load(['neg_epi_','bp_dr_act_p',num2str(patnum),epi_loctag,'.txt']);
epi_drthresh = sort(-epi_drthresh(:,1));

sub_thdata = load(['neg_sub_','bp_dc_act_p',num2str(patnum),sub_loctag,'.txt']);
sub_thdata_Lbox = zeros(10,10);
sub_thdata_Rbox = zeros(10,10);
sub_drthresh = load(['neg_sub_','bp_dr_act_p',num2str(patnum),sub_loctag,'.txt']);
sub_drthresh = sort(-sub_drthresh(:,1));

cd ..;

cd(['patient',num2str(patnum),'_potentials']);
epi_iunit = load(['epi_bp_current_p',num2str(patnum),epi_loctag,'.txt']);
sub_iunit = load(['sub_bp_current_p',num2str(patnum),sub_loctag,'.txt']);
cd ..;

epi_thdata_Lbox(:) = -epi_thdata(1:(noaxons/2),1);
epi_thdata_Rbox(:) = -epi_thdata((noaxons/2+1):noaxons,1);
sub_thdata_Lbox(:) = -sub_thdata(1:(noaxons/2),1);
sub_thdata_Rbox(:) = -sub_thdata((noaxons/2+1):noaxons,1);

% 1 V/A (ohm)
epi_ra = 1e3;%/epi_iunit; 
sub_ra = 1e3;%/sub_iunit;

xper = 10;
epi_drxper = epi_drthresh(round(noaxons/xper))/epi_ra/1e-3;
sub_drxper = sub_drthresh(round(noaxons/xper))/sub_ra/1e-3;

% ----------------------- Plot Monpolar Case ------------------------------

t11_vertlv = {'T12','L1','L2','L3','L4','L5','S1','S2','S3','S4/5'};
t8_vertlv  = {'T9','T10/11','T12','L1','L2/3','L3/4','L5','S1/2','S2/3','S4/5'};
nolv = length(t8_vertlv);

epi_y1 = 0.8*min(min(epi_thdata_Lbox/epi_ra/1e-3));
epi_y2 = 1.2*max(max(epi_thdata_Lbox/epi_ra/1e-3));
sub_y1 = 0.8*min(min(sub_thdata_Lbox/sub_ra/1e-3));
sub_y2 = 1.2*max(max(sub_thdata_Lbox/sub_ra/1e-3));

x_lv = 0:1:(nolv+1);

figure;
subplot(2,1,1);
hold on;
h1 = boxplot(epi_thdata_Lbox/epi_ra/1e-3,'labels',t8_vertlv,...
        'Color','k');
plot(x_lv,epi_drxper*ones(size(x_lv)),'k--','LineWidth',3);
ttl = ['Epidural Stimulation (Patient ',num2str(patnum),'): ',...
       num2str(th_val),'^o Lateral Deviation'];
title(ttl,'FontSize',30,'FontWeight','b');
%ylim([epi_y1,epi_y2]);
%ylim([5,18]);
set(gca,'FontSize',24);
set(gca,'xtick',1:size(epi_thdata_Lbox,2),'xticklabel',t8_vertlv,'FontSize',30);
set(h1,'LineWidth',3);
legend([num2str(xper),' % DR activation'],'Location','NE');
hold off;

subplot(2,1,2);
hold on;
h2 = boxplot(sub_thdata_Lbox/sub_ra/1e-3,'labels',t8_vertlv,'color','k');
plot(x_lv,sub_drxper*ones(size(x_lv)),'k--','LineWidth',3);
ttl = ['Intradural Stimulation (Patient ',num2str(patnum),'): ',...
       num2str(th_val),'^o Lateral Deviation'];
title(ttl,'FontSize',30,'FontWeight','b');
%ylim([sub_y1,sub_y2]);
%ylim([0.01,6]);
set(gca,'FontSize',24);
set(gca,'xtick',1:size(sub_thdata_Lbox,2),'xticklabel',t8_vertlv,'FontSize',30);
set(h2,'LineWidth',3);
hold off;

set(gcf,'NextPlot','add');
axes;
h = ylabel({'Stimulation Threshold Current (mA)',''},'FontSize',30);
set(gca,'Visible','off');
set(h,'Visible','on');
